#ifndef __BAYESIANCLUSTERING_H__
#define __BAYESIANCLUSTERING_H__

#include <cstdio>
#include <string>
#include <ctime>
#include <cmath>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <iomanip>

#include "alignment.h"

using namespace std;

//Global Variables
//*******************************
#define MAXSize 20000
extern int N;
extern int Step;
extern int k;
extern double* Distrk;
extern double* Time;
//extern int* Indicator;
extern double* Pi;
extern double* PiPrev;
extern int* Center;
extern double* Sigma;
extern bool* SigmaChange;
extern double* Likelihood;
extern double* LikelihoodRatio; //for UpdateLikelihood;
//************************************
extern double** LikelihoodTemp;
//************************************
extern char** SeqID;     //for Sequence ID
extern char** seq;  //for reading Fasta
extern char*** EigenSeqs; //for further iteration use
extern float** X;         //Distance Matrix
//double X2[MAXSize][MAXSize];
extern int* Z;       //missing data
extern int* NN;       //NN[i-1]=# of samples in i-th cluster
extern int* IterNN;   //Record # of seqs in a cluster in this round (everyone is a center seq in previous round)
extern int* Abundance; //for further iteration use
extern double* Centerlikelihood;  //square sum of current center
extern gsl_rng *r;
extern bool* Weight;  //for Birth-Death control;
extern unsigned long int WeightSum;
//*******************************

//Result Process
//*******************************
extern int Resultk;
extern double* ResultPi;
extern int* ResultCenter;
extern double* ResultSigma;
extern double* ResultLikelihood;
//*******************************

/*
//Golden Standard Store
//*******************************
extern int StandardZ[MAXSize];
extern int Standardk;
extern int StandardNN[MAXSize];
extern double Purity,Pu;
extern double VI,V;
extern int EditDistance;
extern double ED;
*/

//Hyperparameters
//*******************************
//For Pi
extern double* Gamma;
//For Sigma
extern const double Alpha;
extern double Beta;
extern const double g;
extern double h;
extern double SigmaSum;
extern double Lower, Upper;
//For New-born Pi
extern double Gamma1;
//For Birth-rate
extern double Lamdab;
//For Stationary Time
extern double Lamda;
//*******************************
//Result Storage
//********************************
typedef struct _BayesianClusteringResult{
        int ClusterNumber;
        char** ClusterCenterID;
        char** ClusterCenterSeq;
        double* ClusterStandardDeviation;
        int* ClusterSize;
        int* ClusterIterationSize;
        string** ClusterMember;
        string** ClusterMemberSeq;
        bool** ClusterMemberFlag;        
}BayesianClusteringResult;

//Functions
//*******************************
void InitRandomGenerator();
void LoadDistanceMatrix(const char*, int);
int GenerateInitial(int);
int UpdateZ();
int UpdatePi();
int UpdateCenter(int, bool);
int UpdateSigma();
int UpdateLikelihoodRatio();
int Birth(int);
int Death(int,int);
int ProcessSampling(int, int);
int ResultProcess();
//int CriteriaCalc(int);
double factorial(int,int);
double gaussian(double, double);
BayesianClusteringResult BayesianClustering(const char*,int, bool, BayesianClusteringResult&, int);
void NewMem(int);
void DeleteMem(int,int);
char** readBlockSeq(char*, unsigned int*, int);
extern "C"
{
    char** readFastASeq(char*, unsigned int*);
    #define DELETE(a) \
    {\
        while((a)!=NULL)\
        {\
		    delete[] (a);\
		    (a) = NULL;\
	    }\
    }
}
//*******************************

//Time Inspection
//*******************************
extern double TimeZ;
extern double TimeCenter;
extern double TimeSigma;
extern double TimePi;
extern double TimeBirth;
extern double TimeDeath;
extern double TimeLikelihood;



#endif//__BAYESIANCLUSTERING_H__
